{
 "metadata": {
  "name": "",
  "signature": "sha256:aa135613c0bb1b4755ec892c741d6881277cd5f5a3287d533276219142c31188"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Variational Inference for Gaussian Process Classification Models"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "By Wu Lin - <a href=\"mailto:yorker.lin@gmail.com\">yorker.lin@gmail.com</a> - <a href=\"https://github.com/yorkerlin\">https://github.com/yorkerlin</a> Suggested by Heiko Strathmann and Mohammad Emtiyaz Khan\n",
      "\n",
      "Based on the notebook of Gaussian Processes by Heiko Strathmann - <a href=\"mailto:heiko.strathmann@gmail.com\">heiko.strathmann@gmail.com</a> - <a href=\"https://github.com/karlnapf\">github.com/karlnapf</a> - <a href=\"http://herrstrathmann.de\">herrstrathmann.de</a>, and the GP framework of the Google summer of code 2014 project of Wu Lin,  - Google summer of code 2013 projec of Roman Votyakov - <a href=\"mailto:votjakovr@gmail.com\">votjakovr@gmail.com</a> - <a href=\"https://github.com/votjakovr\">github.com/votjakovr</a>, and the Google summer of code 2012 project of Jacob Walker - <a href=\"mailto:walke434@gmail.com\">walke434@gmail.com</a> - <a href=\"https://github.com/puffin444\">github.com/puffin444</a> "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook describes how to do <a href=\"http://en.wikipedia.org/wiki/Variational_Bayesian_methods\">variational inference</a> for <a href=\"http://en.wikipedia.org/wiki/Gaussian_process\">Gaussian Process <a href=\"http://en.wikipedia.org/wiki/Statistical_classification\">classification</a> (GPC) models in Shogun. In a GPC model, a set of GP functions are used to predict the label of a data point given its features. \n",
      "We assume that readers have some background in <a href=\"http://en.wikipedia.org/wiki/Bayesian_statistics\">Bayesian statistics</a>. For the background, please see the <a href=\"http://www.shogun-toolbox.org/static/notebook/current/gaussian_processes.html\">notebook</a> about Gaussian Process. After providing a semi-formal introduction, we illustrate how to do training, make predictions, and automatically learn hyper-parameters for GPC models in Shogun"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "# import all shogun classes\n",
      "from modshogun import *\n",
      "\n",
      "# import all required libraries\n",
      "import scipy\n",
      "import scipy.io\n",
      "import numpy as np\n",
      "from math import exp,sqrt,log\n",
      "import random\n",
      "import time\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Brief Introduction"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In binary classification, we get binary labels in form of a column vector, $\\mathbf{y}\\in\\mathcal{Y}^n=\\{-1,+1\\}^n$, and features as a matrix, $\\mathbf{X}\\in\\mathcal{R}^{n\\times d}$, where $n$ is the number of data points and $d$ is the number of features.  A likelihood function $p(\\mathbf{y}|\\text{f},\\mathbf{X})=p(\\mathbf{y}|\\mathbf{f})$ is used to model data with labels, where a function $\\text{f}:\\mathcal{R}^{d}\\to \\mathcal{R} $ is drawn from a Gaussian Process prior and $\\text{f}(\\mathbf{X})=\\mathbf{f} \\in \\mathcal{R}^{n}$ is a column vector by applying f to each row of $\\mathbf{X}$. Note that the difference between f and $\\mathbf{f}$ is that f is a function while $\\mathbf{f}$ is an n-dimensional vector. $p(\\text{f})$ and $p(\\mathbf{f})$ are also different because $p(\\text{f})$ is a distribution in an infinite dimensional (function) space while $p(\\mathbf{f})$ is a distribution in an finite dimensional space.\n",
      "\n",
      "Given data with labels, our goal is to train a Gaussian Process classifier to fit the data well.\n",
      "In other words, we want to learn posterior, $p(\\text{f}|\\mathbf{y},\\mathbf{X}) \\propto p(\\mathbf{y}| \\text{f},\\mathbf{X}) \\cdot p(\\text{f}) $, given training data points. Note that $p(\\text{f})$ is the prior over the GP functions. In fact, we will see in the next session that all we need is to learn $p(\\mathbf{f}|\\mathbf{y})$ (Note that we use $\\mathbf{f}$ as $\\text{f}(\\mathbf{X})$ for short).\n",
      "\n",
      "The key intuition of variational inference is to approximate the distribution of interest (here $p(\\mathbf{f}|\\mathbf{y})$, which is a non-Gaussian distribution) with a more tractable distribution (here a Gaussian $q(\\mathbf{f}|\\mathbf{y}))$, via minimizing the <a href=\"http://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence\">Kullback\u2013Leibler divergence</a> (KL divergence) \n",
      "\n",
      "$${\\mathrm{KL}}(Q\\|P) = \\int_{-\\infty}^\\infty \\ln\\left(\\frac{q(\\mathbf{f}|\\mathbf{y})}{p(\\mathbf{f}|\\mathbf{y})}\\right) q(\\mathbf{f}|\\mathbf{y}) \\ {\\rm d}\\mathbf{f}$$\n",
      "\n",
      "between the true distribution and the approximation.\n",
      "\n",
      "Please see the next session if you want to know the reason why non-Gaussian $p(\\mathbf{f}|\\mathbf{y})$ cannot be used in inference process of GPC models.\n",
      "\n",
      "Throughout this notebook, we deal with binary classification using <a href=\"http://en.wikipedia.org/wiki/Logit\">inverse-logit</a>, (also known as Bernoulli-logistic function)\n",
      "likelihood to model labels, given by \n",
      "\n",
      "$p(\\mathbf{y}|\\mathbf{f})=\\prod_{i=1}^n p(y_\\text{i}|\\text{f}(\\mathbf{X}_\\text{i}))=\\prod_{i=1}^n \\frac{1}{1-\\exp(-y_\\text{i} \\mathbf{f}_\\text{i})}$,\n",
      "where $\\mathbf{f}_\\text{i}$ is the i-th row of $\\mathbf{X}$ and $\\text{f}(\\mathbf{X}_\\text{i})=\\mathbf{f}_\\text{i}$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "More Detailed Background about GPC Models (Skip if you just want code examples)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In GPC models, our goal is to predict the label ($y_\\text{new}$) of a new data point, a column vector, $\\mathbf{x}_\\text{new}\\in\\mathcal{R}^{d}$ based on $p(y_\\text{new}|\\mathbf{y},\\mathbf{X},\\mathbf{x}_\\text{new})$. \n",
      "According to <a href=\"http://en.wikipedia.org/wiki/Bayes%27_theorem\">Bayes' theorem</a>, we know that <a href=\"http://en.wikipedia.org/wiki/Posterior_predictive_distribution\"> the Posterior predictive distribution</a> is    $$p(y_\\text{new}|\\mathbf{y},\\mathbf{X},\\mathbf{x}_\\text{new})= \\int {p(y_\\text{new}|f,\\mathbf{x}_\\text{new})p(f|\\mathbf{y},\\mathbf{X}) df}$$ \n",
      "Informally, according to <a href=\"http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Marginal_distributions\">the property</a> of GP about marginalization, we have:\n",
      "\n",
      "$$\\int {p(\\mathbf{y}_\\text{new}|\\text{f},\\mathbf{x}_\\text{new})p(\\text{f}|\\mathbf{y},\\mathbf{X}) d\\text{f}}= \\int {p(\\mathbf{y}_\\text{new}|\\mathbf{f}_\\text{new})p(\\mathbf{f}_\\text{new}|\\mathbf{y},\\mathbf{X}) d\\mathbf{f}_\\text{new}}$$.\n",
      "\n",
      "where $\\text{f}(\\mathbf{x}_\\text{new})=\\mathbf{f}_\\text{new}$, $p(\\mathbf{y}_\\text{new}|\\mathbf{f}_\\text{new})=p(\\mathbf{y}_\\text{new}|\\text{f},\\mathbf{x}_\\text{new})$ and $p(\\mathbf{f}_\\text{new}|\\mathbf{y},\\mathbf{X})=p(\\mathbf{f}_\\text{new}|\\mathbf{y},\\mathbf{X}, \\mathbf{x}_\\text{new})$ .\n",
      "\n",
      "The key difference here is that $p(f|\\mathbf{y}, \\mathbf{X})$ is a distribution in infinite dimensional space while $p(\\mathbf{f}_{new}|\\mathbf{y},\\mathbf{X})$ is a distribution in one-dimensional space.\n",
      "\n",
      "Note that $p(\\mathbf{f}_\\text{new}|\\mathbf{y},\\mathbf{X})=\\int {p(\\mathbf{f}_\\text{new}|\\text{f}) p(\\text{f}|\\mathbf{y},\\mathbf{X}) d\\text{f}}$.\n",
      "Similarly, $\\int {p(\\mathbf{f}_\\text{new}|\\text{f}) p(\\text{f}|\\mathbf{y},\\mathbf{X}) d\\text{f}}=\\int {p(\\mathbf{f}_\\text{new}|\\mathbf{f}) p(\\mathbf{f}|\\mathbf{y}) d\\mathbf{f}}$.\n",
      "Again, $p(\\mathbf{f}|\\mathbf{y})=p(\\mathbf{f}|\\mathbf{y}, \\mathbf{X})$ is a distribution in n-dimensional space while $p(\\text{f}|\\mathbf{y},\\mathbf{X})$ is distribution in infinite dimensional space.\n",
      "Informally, according to the definition of GP, the following holds:\n",
      "\n",
      "$p(\\mathbf{f}_\\text{new}|\\mathbf{f})=\\frac{p(\\text{f}(\\mathbf{[\\mathbf{X};\\mathbf{x}_\\text{new}^T]}))}{p(\\text{f}(\\mathbf{\\mathbf{X}}))}$, which follows from the conditional properties of jointly Gaussian variables, where $[\\mathbf{X};\\mathbf{x}_\\text{new}^T]\\in \\mathcal{R}^{(n+1)\\times d}$ and $\\text{f}([\\mathbf{X};\\mathbf{x}_\\text{new}^T]) \\in \\mathcal{R}^{n+1} $\n",
      "\n",
      "Note that If $p(\\mathbf{f}|\\mathbf{y})$ is a Gaussian distribution, $p(\\mathbf{f}_{new}|\\mathbf{y},\\mathbf{X})$ has a close form. However,in classification $p(\\mathbf{f}|\\mathbf{y})$ usually is NOT a Gaussian distribution since $p(\\mathbf{f}|\\mathbf{y}) \\propto p(\\mathbf{y}|\\mathbf{f})p(\\mathbf{f})$, where $p(\\mathbf{f})$ is a Gaussian distribution but $p(\\mathbf{y}|\\mathbf{f}))$ (the likelihood) is a non-Gaussian distribution.\n",
      "\n",
      "Variational inference in GPC is to approximate $p(\\mathbf{f}|\\mathbf{y})$ using a Gaussian distribution, $q(\\mathbf{f}|\\mathbf{y})$, via minimizing the KL divergence, ${\\mathrm{KL}}(Q\\|P)$. Reader may note that KL divergence is asymmetric. If we minimize ${\\mathrm{KL}}(P\\|Q)$ instead of ${\\mathrm{KL}}(Q\\|P)$, it is called <a href=\"http://en.wikipedia.org/wiki/Expectation_propagation\">expectation propagation</a> (EP). "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Variational Inference in GPC Models (Skip if you just want code examples)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As mentioned in the first section, variaitonal inference in GPC models is about **minimizing** the KL divergence given as below:\n",
      "\n",
      "${\\mathrm{KL}}(Q\\|P) = \\int_{-\\infty}^\\infty \\ln\\left(\\frac{q(\\mathbf{f}|\\mathbf{y})}{p(\\mathbf{f}|y)}\\right) q(\\mathbf{f}|\\mathbf{y}) \\ {\\rm d}\\mathbf{f} = \\int_{-\\infty}^\\infty \\ln\\left(\\frac{q(\\mathbf{f}|\\mathbf{y})}{  \\frac{p(\\mathbf{y}|\\mathbf{f})p(\\mathbf{f})}{p(\\mathbf{y})}  }\\right) q(\\mathbf{f}|\\mathbf{y}) \\ {\\rm d}\\mathbf{f}=\\int_{-\\infty}^\\infty \\ln\\left(\\frac{q(\\mathbf{f}|\\mathbf{y})}{ p(\\mathbf{y}|\\mathbf{f})p(\\mathbf{f}|\\mathbf{y})  }\\right) q(\\mathbf{f}|\\mathbf{y}) \\ {\\rm d}\\mathbf{f} - \\text{const}$. \n",
      "\n",
      "\n",
      "Another way to explain variational inference in GPC is to **maximize** a lower bound of the log of marginal likelihood, $\\ln (p(\\mathbf{y}|\\mathbf{X})) $.\n",
      "\n",
      "$\\ln (p(\\mathbf{y}|\\mathbf{X})) = \\ln (\\int_{-\\infty}^\\infty {p(\\mathbf{y}|\\text{f})p(\\text{f}|\\mathbf{X})} d\\text{f}) =  \\ln (\\int_{-\\infty}^\\infty {p(\\mathbf{y}|\\mathbf{f})p(\\mathbf{f})} d\\mathbf{f})= \\ln (\\int_{-\\infty}^\\infty { q(\\mathbf{f}) \\frac{ p(\\mathbf{y}|\\mathbf{f})p(\\mathbf{f})} {q(\\mathbf{f}|\\mathbf{y})}} d\\mathbf{f}) \\geq  (\\int_{-\\infty}^\\infty { q(\\mathbf{f}|\\mathbf{y}) \\ln ( \\frac{ p(\\mathbf{y}|\\mathbf{f})p(\\mathbf{f})} {q(\\mathbf{f}|\\mathbf{y})}} ) d\\mathbf{f})$\n",
      "\n",
      "\n",
      "where the inequality is based on <a href=\"http://en.wikipedia.org/wiki/Jensen%27s_inequality\">Jensen\u2019s inequality</a>.\n",
      "\n",
      "$\\int_{-\\infty}^\\infty { q(\\mathbf{f}|\\mathbf{y}) \\ln ( \\frac{ p(\\mathbf{y}|\\mathbf{f})p(\\mathbf{f})} {q(\\mathbf{f}|\\mathbf{y})}} ) d\\mathbf{f}=- \\int_{-\\infty}^\\infty \\ln\\left(\\frac{q(\\mathbf{f}|\\mathbf{y})}{ p(\\mathbf{y}|\\mathbf{f})p(\\mathbf{f})  }\\right) q(\\mathbf{f}|\\mathbf{y}) \\ {\\rm d}\\mathbf{f} = -E_q[\\ln(q(\\mathbf{f}|\\mathbf{y}))] + E_q[\\ln(p(\\mathbf{f}))] + E_q[\\ln(p(\\mathbf{y}|\\mathbf{f}))]$, \n",
      "\n",
      "where $E_q(\\cdot)$ is the expectation with respect to $q(\\mathbf{f})$.\n",
      "\n",
      "Note that the last term, $E_q[\\ln(p(\\mathbf{y}|\\mathbf{f}))]$, in GPC usually does NOT have a close form.\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Dealing with the Non-closed Form Issue in GPC Models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<table>\n",
      "<tr>\n",
      "<th>Approach to deal with non-closed Form</th>\n",
      "<th>Implemenation in Shogun</th>\n",
      "<th>Supported Inference Method</th>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"http://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature\">Gaussian Numerical Integration for One-dimensional Space</a></td>\n",
      "<td><a href=\"http://www.shogun-toolbox.org/doc/en/3.0.0/classshogun_1_1CIntegration.html\">Integration</a>, <a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLogitVGLikelihood.html\">LogitVGLikelihood (Bernoulli-logistic function)</a> </td>\n",
      "\n",
      "\n",
      "<td>\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSingleLaplacianInferenceMethodWithLBFGS.html\">Laplace Method With LBFGS</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCovarianceInferenceMethod.html\">Covariance Variational Method</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLApproxDiagonalInferenceMethod.html\">Mean-field Variational Method</a>,\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCholeskyInferenceMethod.html\">Cholesky Variational Method</a>\n",
      "</td>\n",
      "\n",
      "\n",
      "\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"http://www.icml-2011.org/papers/376_icmlpaper.pdf\">The Piecewise Variational Bound for Bernoulli-logistic function</a></td>\n",
      "<td><a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLogitVGPiecewiseBoundLikelihood.html\">LogitVGPiecewiseBoundLikelihood</a></td>\n",
      "<td>\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSingleLaplacianInferenceMethodWithLBFGS.html\">Laplace Method With LBFGS</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCovarianceInferenceMethod.html\">Covariance Variational Method</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLApproxDiagonalInferenceMethod.html\">Mean-field Variational Method</a>,\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCholeskyInferenceMethod.html\">Cholesky Variational Method</a>\n",
      "</td>\n",
      "</tr>\n",
      "\n",
      "<tr>\n",
      "<td><a href=\"http://www.cs.cmu.edu/~lafferty/pub/ctm.pdf\">The Bound for Bernoulli-logistic function from Blei's work</a></td>\n",
      "<td><a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLogitDVGLikelihood.html\">LogitDVGLikelihood</a></td>\n",
      "\n",
      "\n",
      "<td>\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLDualInferenceMethod.html\">Dual Variational Method</a>\n",
      "</td>\n",
      "</tr>\n",
      "\n",
      "\n",
      "<tr>\n",
      "<td><a href=\"https://circle.ubc.ca/handle/2429/43640\">The Jaakkola Bound for Bernoulli-logistic function</a></td>\n",
      "<td>NA</td>\n",
      "<td>NA</td>\n",
      "</tr>\n",
      "\n",
      "<tr>\n",
      "<td><a href=\"https://circle.ubc.ca/handle/2429/43640\">The Bohning Bound for Bernoulli-logistic function</a></td>\n",
      "<td>NA</td>\n",
      "<td>NA</td>\n",
      "</tr>\n",
      "\n",
      "\n",
      "<tr>\n",
      "<td><a href=\"http://arxiv.org/abs/1206.6430\">Monte Carlo Method</a></td>\n",
      "<td>NA</td>\n",
      "<td>NA</td>\n",
      "</tr>\n",
      "\n",
      "</table>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can show that $E_q[\\ln(p(\\mathbf{y}|\\mathbf{f}))]$ can be expressed in term of summation of one-dimensional integrations via the similar <a href=\"http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Marginal_distributions\">marginalization property</a> of multivariate Gaussian distribution. Note that according to <a href=\"http://en.wikipedia.org/wiki/De_Finetti%27s_theorem\"> De Finetti's theorem</a> labels ($\\mathbf{y}$) are <a href=\"http://en.wikipedia.org/wiki/Exchangeable_random_variables\"> exchangeable</a> since we assume that these labels are independently generated from latent GP functions.  \n",
      "\n",
      "$$E_q[\\ln(p(\\mathbf{y}|\\mathbf{f}))]=E_q[\\sum_{i=1}^n {\\ln(p(\\mathbf{y}_\\text{i}|\\mathbf{f}_\\text{i})}]=\\sum_{i=1}^n {E_q[\\ln(p(\\mathbf{y}_\\text{i}|\\mathbf{f}_\\text{i})]}=\\sum_{i=1}^n {E_{q_\\text{i}}[\\ln(p(\\mathbf{y}_\\text{i}|\\mathbf{f}_\\text{i})]}$$\n",
      "where $q$ denotes $q(\\mathbf{f}|\\mathbf{y})$ is a multivariable $N(\\mu, \\Sigma)$, $q_i$ denotes $q_i(\\mathbf{f}_\\text{i}|\\mathbf{y}_\\text{i})$ is a univariate $N(\\mu_\\text{i}, \\Sigma_\\text{i,i})$, and $\\mu_\\text{i}$ and $\\Sigma_\\text{i,i}$ are the i-th element of the mean $\\mu$ and the i-th diagonal element of the covariance matrix $\\Sigma$ of $q(\\mathbf{f}|\\mathbf{y})$ respectively.\n",
      "\n",
      "Of course, even $E_{q_\\text{i}}[\\ln(p(\\mathbf{y}_\\text{i}|\\mathbf{f}_\\text{i})]$ usually does not have a close form in GPC. However, it is relatively easy to deal with this one-dimensional problem.  \n",
      "One way to deal with this issue is one dimensional numerical integration.\n",
      "Another way is to use some variational bound. For now, we assume this expectation is given and we will discuss later in the notebook."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Summary of Inference Methods in GPC Models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<table>\n",
      "<tr>\n",
      "<th>Inference Method in GPC</th>\n",
      "<th>Idea of Approximation</th>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSingleLaplacianInferenceMethod.html\">Laplace Method</a></td>\n",
      "<td>Using the second-order Taylor expansion in the mode of the true posterior</td>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCovarianceInferenceMethod.html\">Covariance Variational Method</a></td>\n",
      "<td>Minimizing KL(approximated distribution || true posterior), \n",
      "where the approximated distribution has a complete covariance matrix\n",
      "</td>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLApproxDiagonalInferenceMethod.html\">Mean-field Variational Method</a></td>\n",
      "<td>Minimizing KL(approximated distribution || true posterior), \n",
      "where the approximated distribution has a diagonal covariance matrix</td>\n",
      "</tr>\n",
      "\n",
      "<tr>\n",
      "<td><a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCholeskyInferenceMethod.html\">Cholesky Variational Method</a></td>\n",
      "<td>Minimizing KL(approximated distribution || true posterior), \n",
      "where the approximated distribution has a complete covariance matrix in \n",
      "Cholesky representation</td>\n",
      "</tr>\n",
      "\n",
      "<tr>\n",
      "<td><a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLDualInferenceMethod.html\">Dual Variational Method</a></td>\n",
      "<td>Minimizing KL(approximated distribution || true posterior) in a dual form\n",
      "</td>\n",
      "</tr>\n",
      "\n",
      "<tr>\n",
      "<td><a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CEPInferenceMethod.html\">Expectation Propagation Method</a></td>\n",
      "<td>Minimizing KL(true posterior || approximated distribution)\n",
      "</td>\n",
      "</tr>\n",
      "\n",
      "</table>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Summary of Optimizers used in Inference Methods in GPC Models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<table>\n",
      "<tr>\n",
      "<th>Optimizer</th>\n",
      "<th>Inference Method</th>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"http://en.wikipedia.org/wiki/Limited-memory_BFGS\">Limited-memory BFGS</a></td>\n",
      "<td>\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSingleLaplacianInferenceMethodWithLBFGS.html\">Laplace Method With LBFGS</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCovarianceInferenceMethod.html\">Covariance Variational Method</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLApproxDiagonalInferenceMethod.html\">Mean-field Variational Method</a>,\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLDualInferenceMethod.html\">Dual Variational Method</a>,\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCholeskyInferenceMethod.html\">Cholesky Variational Method</a>\n",
      "</td>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization\">Newton-Raphson</a></td>\n",
      "<td><a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSingleLaplacianInferenceMethod.html\">Laplace Method</a>\n",
      "</td>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"http://en.wikipedia.org/wiki/Conjugate_gradient_method\">Conjugate Gradient</a></td>\n",
      "<td>NA</td>\n",
      "</tr>\n",
      "\n",
      "\n",
      "</table>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Summary of Supported Line Search Methods used in Inference Methods for GPC models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<table>\n",
      "<tr>\n",
      "<th>Line Search Method in Limited-memory BFGS Optimizer</th>\n",
      "<th>Inference Method</th>\n",
      "<th>Value in <a href=\"http://www.shogun-toolbox.org/doc/en/latest/lbfgs_8cpp.html\">the ShogunAPI</a></th>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.30.660\">MoreThuente method</a></td>\n",
      "<td>\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSingleLaplacianInferenceMethodWithLBFGS.html\">Laplace Method With LBFGS</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCovarianceInferenceMethod.html\">Covariance Variational Method</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLApproxDiagonalInferenceMethod.html\">Mean-field Variational Method</a>,\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCholeskyInferenceMethod.html\">Cholesky Variational Method</a>,\n",
      "(For now, **NOT** support <a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLDualInferenceMethod.html\">Dual Variational Method</a>)\n",
      "</td>\n",
      "<td>0</td>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"http://en.wikipedia.org/wiki/Wolfe_conditions\">Backtracking method with the Armijo condition</a></td>\n",
      "\n",
      "\n",
      "<td>\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSingleLaplacianInferenceMethodWithLBFGS.html\">Laplace Method With LBFGS</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCovarianceInferenceMethod.html\">Covariance Variational Method</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLApproxDiagonalInferenceMethod.html\">Mean-field Variational Method</a>,\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLDualInferenceMethod.html\">Dual Variational Method</a>,\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCholeskyInferenceMethod.html\">Cholesky Variational Method</a>\n",
      "</td>\n",
      "<td>1</td>\n",
      "</tr>\n",
      "<tr>\n",
      "<td><a href=\"http://en.wikipedia.org/wiki/Wolfe_conditions\">Backtracking method with the (regular or strong) Wolfe condition</a></td>\n",
      "\n",
      "<td>\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CSingleLaplacianInferenceMethodWithLBFGS.html\">Laplace Method With LBFGS</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCovarianceInferenceMethod.html\">Covariance Variational Method</a>,\n",
      "<a href=\"www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLApproxDiagonalInferenceMethod.html\">Mean-field Variational Method</a>,\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLDualInferenceMethod.html\">Dual Variational Method</a>,\n",
      "<a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CKLCholeskyInferenceMethod.html\">Cholesky Variational Method</a>\n",
      "</td>\n",
      "<td>2 (regular),3 (strong, default line search method)</td>\n",
      "\n",
      "</tr>\n",
      "\n",
      "\n",
      "</table>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "A Toy Example of Variational Inference"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we present a 2-D example for GP classification using variational inference. The gaussian prior and the likelihood are shown in Fig. (a) and (b). For this simple example, we can numerically compute the unnormalized true posterior, shown in Fig. (c), using the Bayes rule: $p(\\mathbf{f}|\\mathbf{y}) \\propto p(\\mathbf{y}|\\mathbf{f})p(\\mathbf{f})$, where $p(\\mathbf{f})$ is the prior and $p(\\mathbf{y}|\\mathbf{f}))$ is the likelihood. The approximated posterior obtained using various methods are shown in Fig. (g)-(i).\n",
      "These figures can be obtained by using the following code. Note that this example is closely related to the example at page 7-8 of <a href=\"http://www.jmlr.org/papers/volume9/nickisch08a/nickisch08a.pdf\">this paper</a>."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Gaussian_points(Sigma,mu,xmin,xmax,ymin,ymax,delta=0.01):\n",
      "    \"\"\"\n",
      "    This function is used to evaluate the likelihood of a Gaussian distribution on mesh.\n",
      "    \n",
      "    Parameters:\n",
      "         Sigma - covariance of Gaussian\n",
      "         mu - mean of Gaussian\n",
      "         xmin, xmax, ymin, ymax, delta - defines mesh\n",
      "         \n",
      "    Returns:\n",
      "    X, Y, Z, where Z = log(p(X,Y)) and p is a bivariate gaussian (mu, Sigma) evaluated at a mesh X,Y\n",
      "    \"\"\"\n",
      "\n",
      "    xlist = np.arange(xmin, xmax, delta)\n",
      "    ylist = np.arange(ymin, ymax, delta)\n",
      "    X, Y = np.meshgrid(xlist, ylist)\n",
      "    model = GaussianDistribution(mu, Sigma)\n",
      "    Z = np.zeros(len(X)*len(Y))\n",
      "    idx = 0\n",
      "    for items in zip(X,Y):\n",
      "        for sample in zip(items[0],items[1]):\n",
      "            sample = np.asarray(sample)\n",
      "            Z[idx] = model.log_pdf(sample)\n",
      "            idx += 1\n",
      "    Z = np.asarray(Z).reshape(len(X),len(Y))\n",
      "    return (X,Y,Z)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def likelihood_points(X,Y,labels,likelihood):\n",
      "    \"\"\"\n",
      "    This function is used to evaluate a given likelihood function on a given mesh of points\n",
      "    \n",
      "    Parameters:\n",
      "        X, Y - coordiates of the mesh\n",
      "        labels - labels for the likelihood\n",
      "        likelihood - likelihood function\n",
      "        \n",
      "    Returns:\n",
      "        Z - log(p(X,Y,labels)), where p comes from likelihood\n",
      "    \"\"\"\n",
      "    \n",
      "    Z = np.zeros(len(X)*len(Y))    \n",
      "    idx = 0\n",
      "    for items in zip(X,Y):\n",
      "        for sample in zip(items[0],items[1]):\n",
      "            sample = np.asarray(sample)\n",
      "            lpdf = likelihood.get_log_probability_f(labels, sample).sum()\n",
      "            Z[idx] = lpdf\n",
      "            idx += 1\n",
      "    Z = np.asarray(Z).reshape(len(X),len(Y))\n",
      "    return Z"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def approx_posterior_plot(methods, kernel_func, features, mean_func, \n",
      "                          labels, likelihoods, kernel_log_scale, \n",
      "                          xmin, xmax, ymin, ymax, delta, plots):\n",
      "    \"\"\"\n",
      "    This function is used to generate points drawn from approximated posterior and plot them\n",
      "        \n",
      "    Parameters:\n",
      "        methods - a list of methods used to approximate the posterior\n",
      "        kernel_func - a covariance function for GP\n",
      "        features - X\n",
      "        mean_func -  mean function for GP\n",
      "        labels - Y\n",
      "        likelihood- a data likelihood to model labels\n",
      "        kernel_log_scale - a hyper-parameter of covariance function\n",
      "        xmin, ymin, xmax, ymax, delta - used in linespace function: maximum and minimum values used to generate points from an approximated posterior\n",
      "        plots\n",
      "    \n",
      "    Returns:\n",
      "        Nothing\n",
      "    \"\"\"\n",
      "    \n",
      "    (rows, cols) = plots.shape\n",
      "    methods = np.asarray(methods).reshape(rows, cols)\n",
      "    likelihoods = np.asarray(likelihoods).reshape(rows, cols)\n",
      "    for r in xrange(rows):\n",
      "        for c in xrange(cols):\n",
      "            inference = methods[r][c]\n",
      "            likelihood = likelihoods[r][c]\n",
      "            inf = inference(kernel_func, features, mean_func, labels, likelihood())\n",
      "            inf.set_scale(exp(kernel_log_scale))\n",
      "            #get the approximated Gaussian distribution\n",
      "            mu = inf.get_posterior_mean()\n",
      "            Sigma = inf.get_posterior_covariance()\n",
      "            #normalized approximated posterior\n",
      "            (X,Y,Z) = Gaussian_points(Sigma, mu, xmin, xmax, ymin, ymax, delta)\n",
      "            plots[r][c].contour(X, Y, np.exp(Z))\n",
      "            plots[r][c].set_title('posterior via %s'%inf.get_name())\n",
      "            plots[r][c].axis('equal')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#a toy 2D example (data)\n",
      "x=np.asarray([sqrt(2),-sqrt(2)]).reshape(1,2)\n",
      "y=np.asarray([1,-1])\n",
      "\n",
      "features = RealFeatures(x)\n",
      "labels = BinaryLabels(y)\n",
      "kernel_log_sigma = 1.0\n",
      "kernel_log_scale = 1.5\n",
      "\n",
      "#a mean function and a covariance function for GP\n",
      "mean_func = ConstMean()\n",
      "#using log_sigma as a hyper-parameter of GP instead of sigma\n",
      "kernel_sigma = 2*exp(2*kernel_log_sigma)\n",
      "kernel_func = GaussianKernel(10, kernel_sigma)\n",
      "kernel_func.init(features, features)\n",
      "#a prior distribution derived from GP via applying the mean function and the covariance function to data\n",
      "Sigma = kernel_func.get_kernel_matrix()\n",
      "Sigma = Sigma * exp(2.0*kernel_log_scale)\n",
      "mu = mean_func.get_mean_vector(features)\n",
      "\n",
      "delta = 0.1\n",
      "xmin = -4\n",
      "xmax = 6\n",
      "ymin = -6\n",
      "ymax = 4\n",
      "\n",
      "#a prior (Gaussian) derived from GP\n",
      "(X,Y,Z1) = Gaussian_points(Sigma, mu, xmin, xmax, ymin, ymax, delta)\n",
      "\n",
      "col_size = 6\n",
      "f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(col_size*3,col_size))\n",
      "\n",
      "ax1.contour(X, Y, np.exp(Z1))\n",
      "ax1.set_title('Prior')\n",
      "ax1.axis('equal')\n",
      "\n",
      "#likelihood (inverse logit, A.K.A. Bernoulli-logistic)\n",
      "#likelihoods classes for inference methods  (see the table above)\n",
      "likelihoods=[\n",
      "LogitLikelihood,\n",
      "LogitLikelihood,\n",
      "LogitVGLikelihood,\n",
      "LogitVGLikelihood,\n",
      "LogitVGLikelihood,\n",
      "LogitDVGLikelihood\n",
      "]\n",
      "\n",
      "\n",
      "Z2 = likelihood_points(X,Y,labels,LogitLikelihood())\n",
      "ax2.contour(X, Y, np.exp(Z2))\n",
      "ax2.set_title('Likelihood')\n",
      "ax2.axis('equal')\n",
      "\n",
      "#a unnormalized true posterior (non-Gaussian)\n",
      "Z3 = Z1+Z2\n",
      "ax3.contour(X, Y, np.exp(Z3))\n",
      "ax3.set_title('True posterior')\n",
      "ax3.axis('equal')\n",
      "\n",
      "f, plots =plt.subplots(2, 3, figsize=(col_size*3,col_size*2))\n",
      "\n",
      "#Inference methods used to approximate a posterior (see the table below)\n",
      "methods=[\n",
      "SingleLaplacianInferenceMethod,\n",
      "SingleLaplacianInferenceMethodWithLBFGS,\n",
      "KLApproxDiagonalInferenceMethod,\n",
      "KLCovarianceInferenceMethod,\n",
      "KLCholeskyInferenceMethod,\n",
      "KLDualInferenceMethod\n",
      "]\n",
      "\n",
      "approx_posterior_plot(methods, kernel_func, features, mean_func, labels, likelihoods, kernel_log_scale, \n",
      "                      xmin, xmax, ymin, ymax, delta, plots)\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remark:\n",
      "The true posterior is non-Gaussian and the normalized constant usually is difficult to compute.\n",
      "All approximated posteriors are Gaussian. A posterior Gaussian distribution is required to perform inference and predictions.\n",
      "The Laplace approximation is not ideal in term of approximating a posterior and making predictions although it is fast.\n",
      "With the price of speed, variational methods such as the KLCovariance and the KLCholesky offer more accurate approximations compared than the Laplace method. The KLDual method also can have a good approximation and in many cases are better than the Laplace method. On the other hand, the KLApproxDiagonal (mean-field) method is as fast as the Laplace method and for certain datasets, the KLApproxDiagonal method is more accurate than Laplace method in term of making predictions. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "A Real-World Example"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We apply variational methods to the sonar data set. The data set can be found at <a href=\"https://archive.ics.uci.edu/ml/machine-learning-databases/undocumented/connectionist-bench/sonar/\">here</a>. This is a binary classification problem. The goal of this example is to compare all inference methods implemented in shogun in term of making prediction predictions.  The following is the description of the data set. This example is about reproducing the example at page 29 of <a href=\"http://www.jmlr.org/papers/volume9/nickisch08a/nickisch08a.pdf\">this paper</a>. Please refer the paper about the definition of the information bit measure.\n",
      "\n",
      "The class \"sonar.mines\" contains 111 data points  obtained by bouncing sonar\n",
      "signals off a metal cylinder at various angles and under various\n",
      "conditions.  The class \"sonar.rocks\" contains 97 data points obtained from\n",
      "rocks under similar conditions.  The transmitted sonar signal is a\n",
      "frequency-modulated chirp, rising in frequency.  The data set contains\n",
      "signals obtained from a variety of different aspect angles, spanning 90\n",
      "degrees for the cylinder and 180 degrees for the rock.\n",
      "Each data point is a set of 60 features in the range 0.0 to 1.0.  Each feature\n",
      "represents the energy within a particular frequency band, integrated over\n",
      "a certain period of time.  The integration aperture for higher frequencies\n",
      "occur later in time, since these frequencies are transmitted later during\n",
      "the chirp.\n",
      "The label associated with each record contains the letter \"R\" if the object\n",
      "is a rock and \"M\" if it is a mine (metal cylinder).  The features in the\n",
      "labels are in increasing order of aspect angle, but they do not encode the\n",
      "angle directly."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def learning_example1(inference, linesearch, likelihood, x_train, y_train, x_test, y_test, kernel_log_sigma, kernel_log_scale):\n",
      "    \n",
      "    \"\"\"\n",
      "    This function is used to train a GP classifer given an inference method\n",
      "        \n",
      "    Parameters:\n",
      "        inference - an inference methods used to train the classifer\n",
      "        linesearch - a linesearch method used in the inference method \n",
      "        likelihood - likelihood function to model data\n",
      "        x_train, x_test - X for training and testing\n",
      "        y_train, y_test - Y for training and testing\n",
      "        kernel_log_sigma, kernel_log_scale hyper-parameters of covariance function    \n",
      "    Returns:\n",
      "        predictive result of the testing data set\n",
      "    \"\"\"\n",
      "\n",
      "    mean_func = ZeroMean()\n",
      "    kernel_sigma = 2*exp(2*kernel_log_sigma);\n",
      "    kernel_func = GaussianKernel(10, kernel_sigma)\n",
      "\n",
      "    #Y is a sample-by-1 vector\n",
      "    labels_train = BinaryLabels(y_train)\n",
      "    labels_test = BinaryLabels(y_test)\n",
      "    #X is a feature-by-sample matrix\n",
      "    features_train=RealFeatures(x_train)\n",
      "    features_test=RealFeatures(x_test)\n",
      "\n",
      "    inf = inference(kernel_func, features_train, mean_func, labels_train, likelihood)\n",
      "    inf.set_scale(exp(kernel_log_scale))\n",
      "    try:\n",
      "        #setting lbfgs parameters\n",
      "        inf.set_lbfgs_parameters(100,2000,int(linesearch),2000)\n",
      "        #used to make sure the kernel matrix is positive definite\n",
      "        inf.set_noise_factor(1e-6)\n",
      "        inf.set_min_coeff_kernel(1e-5)\n",
      "        inf.set_max_attempt(10)\n",
      "    except:\n",
      "        pass\n",
      "\n",
      "    gp = GaussianProcessClassification(inf)\n",
      "    gp.train()\n",
      "    prob=gp.get_probabilities(features_test)\n",
      "\n",
      "    return prob, inf.get_name()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "labels={\n",
      "    \"m\":1,\n",
      "    \"r\":-1,\n",
      "       }\n",
      "\n",
      "def extract(inf, train_size):\n",
      "    \"\"\"\n",
      "    This function is used to processing raw data\n",
      "        \n",
      "    Parameters:\n",
      "        inf - the path of the raw data\n",
      "        train_size - number of data points used for testing   \n",
      "        \n",
      "    Returns:\n",
      "        x_train, x_test - X for training and testing\n",
      "        y_train, y_test - Y for training and testing\n",
      "        B - for computing the information bit\n",
      "    \"\"\"\n",
      "\n",
      "    random.seed(1)\n",
      "    x=[]\n",
      "    y=[]\n",
      "    for line in open(inf):\n",
      "        line=line.strip()\n",
      "        info=line.split(',')\n",
      "        label=labels[info[-1].lower()]\n",
      "        x.append(map(float, info[:-1]))\n",
      "        y.append(label)\n",
      "        \n",
      "    #train_size should be less than the size of all available data points\n",
      "    assert train_size < len(x) \n",
      "    idx=range(len(y))\n",
      "    random.shuffle(idx)\n",
      "    train_idx=set(idx[:train_size])\n",
      "    test_idx=set(idx[train_size:])\n",
      "    x_train = np.asarray([value for (idx, value) in enumerate(x) if idx in train_idx]).T\n",
      "    y_train = np.asarray([label for (idx, label) in enumerate(y) if idx in train_idx])\n",
      "    x_test = np.asarray([value for (idx, value) in enumerate(x) if idx in test_idx]).T\n",
      "    y_test = np.asarray([label for (idx, label) in enumerate(y) if idx in test_idx])\n",
      "\n",
      "    y_train_positive=(y_train==1).sum()\n",
      "    y_train_negative=(y_train==-1).sum()\n",
      "    y_test_positive=(y_test==1).sum()\n",
      "    y_test_negative=(y_test==-1).sum()\n",
      "\n",
      "    prb_positive=float(y_train_positive)/len(y_train)\n",
      "    B=float(y_test_positive)*log(prb_positive,2)+float(y_test_negative)*log(1.0-prb_positive,2)\n",
      "    B=-B/len(y_test)\n",
      "\n",
      "    return x_train, y_train, x_test, y_test, B"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def approx_bit_plot(methods, linesearchs, likelihoods, x_train, y_train, x_test , y_test, plots, lScale, lSigma, B):\n",
      "    \"\"\"\n",
      "    This function is used to plot information bit figures\n",
      "        \n",
      "    Parameters:\n",
      "        methods - a list of inference methods used to train the classifer\n",
      "        linesearchs - a list of linesearch methods used in the inference method\n",
      "        likelihood - a list of likelihood functions to model data\n",
      "        x_train, x_test - X for training and testing\n",
      "        y_train, y_test - Y for training and testing\n",
      "        plots  - plot objects\n",
      "        lScale, lSigma  - a list of hyper-parameters of covariance function    \n",
      "        B - for computing information bits\n",
      "        \n",
      "    Returns:\n",
      "        predictive result of the testing data set\n",
      "    \"\"\"\n",
      "    #Note that information bit is used to measure effectiveness of an inference method\n",
      "    if len(plots.shape)==1:\n",
      "        rows=1\n",
      "        cols=plots.shape[0]\n",
      "    else:\n",
      "        (rows, cols) = plots.shape\n",
      "\n",
      "    methods=np.asarray(methods).reshape(rows, cols)\n",
      "    linesearchs=np.asarray(linesearchs).reshape(rows, cols)\n",
      "    likelihoods=np.asarray(likelihoods).reshape(rows, cols) \n",
      "\n",
      "    for r in xrange(rows):\n",
      "        for c in xrange(cols):\n",
      "            inference = methods[r][c]\n",
      "            linesearch = linesearchs[r][c]\n",
      "            likelihood = likelihoods[r][c]\n",
      "            try:\n",
      "                likelihood.set_noise_factor(1e-15)\n",
      "                likelihood.set_strict_scale(0.01)\n",
      "            except:\n",
      "                pass\n",
      "            scores=[]\n",
      "            for items in zip(lScale, lSigma):\n",
      "                for parameters in zip(items[0],items[1]):\n",
      "                    lscale=parameters[0]\n",
      "                    lsigma=parameters[1]\n",
      "                    (prob, inf_name)=learning_example1(inference, linesearch, likelihood, x_train, y_train, x_test,  y_test, \n",
      "                                              lscale, lsigma)\n",
      "                    #compute the information bit\n",
      "                    score=0.0\n",
      "                    for (label_idx, prb_idx) in zip(y_test, prob):\n",
      "                        if label_idx==1:\n",
      "                            score+=(1.0+label_idx)*log(prb_idx,2)\n",
      "                        else: #label_idx==-1:\n",
      "                            score+=(1.0-label_idx)*log(1.0-prb_idx,2)\n",
      "                    score=score/(2.0*len(y_test))+B\n",
      "                    scores.append(score)\n",
      "            scores=np.asarray(scores).reshape(len(lScale),len(scores)/len(lScale))\n",
      "\n",
      "            if len(plots.shape)==1:\n",
      "                sub_plot=plots\n",
      "            else:\n",
      "                sub_plot=plots[r]\n",
      "            CS =sub_plot[c].contour(lScale, lSigma, scores)\n",
      "            sub_plot[c].clabel(CS, inline=1, fontsize=10)\n",
      "            sub_plot[c].set_title('information bit for %s'%inf_name)\n",
      "            sub_plot[c].set_xlabel('log_scale')\n",
      "            sub_plot[c].set_ylabel('log_sigma')\n",
      "            sub_plot[c].axis('equal')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#an example for the sonar data set\n",
      "train_size=108\n",
      "(x_train, y_train, x_test, y_test, B)=extract('../../../data/uci/sonar/sonar.all-data', train_size)\n",
      "inference_methods =[\n",
      "                  EPInferenceMethod,\n",
      "                  KLDualInferenceMethod,\n",
      "                  SingleLaplacianInferenceMethod,\n",
      "                  KLApproxDiagonalInferenceMethod,\n",
      "                  KLCholeskyInferenceMethod,\n",
      "                  KLCovarianceInferenceMethod,\n",
      "                  ]\n",
      "\n",
      "likelihoods =[\n",
      "            LogitLikelihood(),\n",
      "            LogitDVGLikelihood(), #KLDual method uses a likelihood class that supports dual variational inference\n",
      "            LogitLikelihood(),\n",
      "            LogitVGLikelihood(), #KL method uses a likelihood class that supports variational inference\n",
      "            LogitVGLikelihood(), #KL method uses a likelihood class that supports variational inference\n",
      "            LogitVGLikelihood(), #KL method uses a likelihood class that supports variational inference\n",
      "            ]\n",
      "\n",
      "linesearchs =[\n",
      "            3,\n",
      "            1, #KLDual method using the type 3 line_search method will cause an error during computing the information bit\n",
      "            3,\n",
      "            3,\n",
      "            3,\n",
      "            3,\n",
      "            ]\n",
      "\n",
      "col_size=8\n",
      "lscale_min=0.0\n",
      "lscale_max=5.0\n",
      "lsigma_min=0.0\n",
      "lsigma_max=5.0\n",
      "delta=0.5\n",
      "scale=5.0\n",
      "\n",
      "lscale_list = np.arange(lscale_min, lscale_max, delta)\n",
      "lsigma_list = np.arange(lsigma_min, lsigma_max, delta*scale)\n",
      "lScale, lSigma = np.meshgrid(lscale_list, lsigma_list)\n",
      "f, plots =plt.subplots(2, 3, figsize=(col_size*3,col_size*2))\n",
      "\n",
      "approx_bit_plot(inference_methods, linesearchs, likelihoods, x_train, y_train,  x_test, y_test, plots, lScale, lSigma, B)\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that the information bit is used to measure the accuracy of classification. \n",
      "If there is a perfect classification, the information bit should be 1. \n",
      "From these figures, we can observe that the Laplace method is worst among these methods in term of information bit. The EP method, the KLCholesky method and the KLCovariance method are equally good. Surprisingly, for this data set the KLApproxDiagonal method performs as good as the EP method in term of information bit. Note that the KLApproxDiagonal method enforces sparsity in non-diagonal entries of the posterior co-variance matrix.  "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Another Real-world Example"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We use the US Postal Service (USPS) database of handwritten digits as another example. The dataset refers to numeric data obtained from the scanning of handwritten digits from envelopes by the U.S. Postal Service. The original scanned digits are binary and of different sizes and orientations; the images here have been deslanted and size normalized, resulting in 16 x 16 grayscale images.\n",
      "\n",
      "We consider to classify the images of digit 3 from images of digit 5 in this example, which means only images of digit 3 and digit 5 are used. Note that this is example is also used in session 3.7.3 of the textbook, <a href=\"http://www.gaussianprocess.org/gpml/\">Gaussian Processes for Machine Learning</a>."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def binary_extract(idx, features, labels):\n",
      "    \"\"\"\n",
      "    This function is used to extract a trainning set and a test set\n",
      "    \n",
      "    Parameters:\n",
      "        idx - a 2-wide list of digits (eg, [0,3] means to extract a sub set of images and labels for digit 3 and digit 0) \n",
      "        features - the whole image set \n",
      "        labels - the whole label set\n",
      "        \n",
      "    Returns:\n",
      "        binary_features - X (images for binary classification)\n",
      "        binary_labels - y (labels for binary classification)\n",
      "    \"\"\"\n",
      "    #binary classification\n",
      "    assert len(idx) == 2\n",
      "    positive_idx = (labels[idx[0],:] == 1)\n",
      "    negative_idx = (labels[idx[1],:] == 1)\n",
      "    binary_idx = (positive_idx | negative_idx)\n",
      "    ds = binary_idx.shape[-1]\n",
      "\n",
      "    bin_labels = np.zeros(ds) \n",
      "    bin_labels[positive_idx] = 1\n",
      "    bin_labels[negative_idx] = -1\n",
      "\n",
      "    binary_features = (features[:,binary_idx])\n",
      "    binary_labels = (bin_labels[binary_idx])\n",
      " \n",
      "    positive_count = bin_labels[positive_idx].shape[-1]\n",
      "    negative_count = bin_labels[negative_idx].shape[-1]\n",
      "    binary_count = binary_labels.shape[-1]\n",
      " \n",
      "    print \"There are %d positive samples and %d negative samples\" %(positive_count, negative_count)\n",
      "    return binary_features, binary_labels"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def learning_example2(inference, likelihood, x_train, x_test, y_train, y_test):\n",
      "    \"\"\"\n",
      "    This function is used to train a GP classifer given an inference method\n",
      "        \n",
      "    Parameters:\n",
      "        inference - an inference methods used to train the classifer\n",
      "        likelihood - likelihood function to model data\n",
      "        x_train, x_test - X for training and testing\n",
      "        y_train, y_test - Y for training and testing \n",
      "        \n",
      "    Returns:\n",
      "        Nothing\n",
      "    \"\"\"\n",
      "        \n",
      "    #train a GP classifer\n",
      "    error_eval = ErrorRateMeasure()\n",
      "    mean_func = ConstMean()\n",
      "    #set hyper-parameters of covariance function\n",
      "    kernel_log_sigma = 1.0\n",
      "    kernel_sigma = 2*exp(2*kernel_log_sigma);\n",
      "    kernel_func = GaussianKernel(10, kernel_sigma)\n",
      "\n",
      "    #sample by 1\n",
      "    labels_train = BinaryLabels(y_train)\n",
      "    labels_test = BinaryLabels(y_test)\n",
      "    #feature by sample\n",
      "    features_train=RealFeatures(x_train)\n",
      "    features_test=RealFeatures(x_test)\n",
      "\n",
      "    kernel_log_scale = 1.0\n",
      "\n",
      "    inf = inference(kernel_func, features_train, mean_func, labels_train, likelihood)\n",
      "    print \"\\nusing %s\"%inf.get_name()\n",
      "    \n",
      "    inf.set_scale(exp(kernel_log_scale))\n",
      "    try:\n",
      "        inf.set_lbfgs_parameters(100,80,3,80)\n",
      "    except:\n",
      "        pass\n",
      "\n",
      "    start = time.time()\n",
      "    gp = GaussianProcessClassification(inf)\n",
      "    gp.train()\n",
      "    end = time.time()\n",
      "    print \"cost %.2f seconds at training\"%(end-start)\n",
      "    nlz=inf.get_negative_log_marginal_likelihood()\n",
      "    print \"the negative_log_marginal_likelihood is %.4f\"%nlz\n",
      "    start = time.time()\n",
      "    #classification on train_data\n",
      "    pred_labels_train = gp.apply_binary(features_train)\n",
      "    #classification on test_data\n",
      "    pred_labels_test = gp.apply_binary(features_test)\n",
      "    end = time.time()    \n",
      "    print \"cost %.2f seconds at prediction\"%(end-start)\n",
      "    \n",
      "    error_train = error_eval.evaluate(pred_labels_train, labels_train)\n",
      "    error_test = error_eval.evaluate(pred_labels_test, labels_test)\n",
      "    \n",
      "    print \"Train error : %.2f Test error: %.2f\\n\" % (error_train, error_test)  "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#an example for the USPS data set\n",
      "inf='../../../data/toy/usps/usps_resampled/usps_resampled.mat'\n",
      "data=scipy.io.loadmat(inf)\n",
      "train_labels=data['train_labels']\n",
      "test_labels=data['test_labels']\n",
      "train_features=data['train_patterns']\n",
      "test_features=data['test_patterns']\n",
      "#using images of digit 3 and digit 5 from the dataset\n",
      "idx=[3,5]\n",
      "#Note that \n",
      "#y_train and y_test are followed the definition in the first session \n",
      "#the transpose of x_train and x_test are followed the definition in the first session\n",
      "print \"Training set statistics\"\n",
      "(x_train, y_train)=binary_extract(idx,train_features, train_labels)\n",
      "print \"Test set statistics\"\n",
      "(x_test, y_test)=binary_extract(idx,test_features, test_labels)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Laplace Method in GPC Models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<a href=\"http://en.wikipedia.org/wiki/Laplace%27s_method\">The Laplace method</a> can be viewed as a special case of variational method.\n",
      "The idea of the Laplace method is to use the second-order <a href=\"http://en.wikipedia.org/wiki/Taylor_series\">Taylor  expansion</a> in the <a href=\"http://en.wikipedia.org/wiki/Mode_%28statistics%29\">mode</a>, $\\mathbf{\\hat{f}}$ of $p(\\mathbf{f}|\\mathbf{y})$ to approximate $p(\\mathbf{f}|\\mathbf{y})$, which is given below:\n",
      "\n",
      "$p(\\mathbf{f}|\\mathbf{y}) \\approx \\hat{p}(\\mathbf{f}|\\mathbf{y})$ followed by $N(\\mathbf{\\hat{f}}, H_{\\mathbf{\\hat{f}}}^{-1})$ ,where $H_{\\mathbf{\\hat{f}}}$ is the <a href=\"http://en.wikipedia.org/wiki/Hessian_matrix\">Hessian matrix</a> in $\\mathbf{\\hat{f}}$\n",
      "\n",
      "Therefore, the KL divergence, ${\\mathrm{KL}}(Q\\|P)$, is approximated by $\\int_{-\\infty}^\\infty \\ln\\left(\\frac{q(\\mathbf{f}|\\mathbf{y})}{\\hat{p}(\\mathbf{f}|y)}\\right) q(\\mathbf{f}|\\mathbf{y}) \\ {\\rm d}\\mathbf{f}$.\n",
      "\n",
      "By **minimizing** $\\int_{-\\infty}^\\infty \\ln\\left(\\frac{q(\\mathbf{f}|\\mathbf{y})}{\\hat{p}(\\mathbf{f}|y)}\\right) q(\\mathbf{f}|\\mathbf{y}) \\ {\\rm d}\\mathbf{f}$\n",
      "we can get $q(\\mathbf{f}|\\mathbf{y})= \\hat{p}(\\mathbf{f}|\\mathbf{y})$.\n",
      "In practice, we can use <a href=\"http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization\">Newton-Raphson</a> optimizer or <a href=\"http://en.wikipedia.org/wiki/Quasi-Newton_method\">Quasi-Newton</a> optimizers(ie, <a href=\"http://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm\">Limited-memory Broyden\u2013Fletcher\u2013Goldfarb\u2013Shanno</a> (LBFGS)) to find $\\mathbf{\\hat{f}}$.\n",
      "Since the likelihood is the inverse-logit, we can show that the objective function, $\\ln(p(\\mathbf{f}|\\mathbf{y}))$, is strictly <a href=\"http://en.wikipedia.org/wiki/Concave\">concave</a>, which means in theory these optimizers can find the same global solution $\\mathbf{\\hat{f}}$.\n",
      "We demonstrate how to apply the Laplace method via Newton-Raphson optimizer and LBFGS optimizer in Shogun to the USPS data set as below."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "inference_methods=[\n",
      "                  SingleLaplacianInferenceMethodWithLBFGS,     #using LBFGS optimizer\n",
      "                  SingleLaplacianInferenceMethod,              #using Newton-Raphson optimizer\n",
      "                  ]\n",
      "likelihood = LogitLikelihood()\n",
      "#using default optimizer and default linesearch method\n",
      "for inference in inference_methods:\n",
      "    learning_example2(inference, likelihood, x_train, x_test, y_train, y_test)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remark:  In practice, the Laplace method is the fastest method among all implemented variational methods in Shogun. However, the second-order Taylor approximation in the mode of marginal likelihood is not always a good approximation, which can be observed in the figures in the previous session for the sonar data set."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Covariance Variational Method in GPC Models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This method is mentioned in <a href=\"http://jmlr.org/papers/volume9/nickisch08a/nickisch08a.pdf\">the paper</a> of Nickisch and Rasmussen in 2008, which is called the KL method in the paper.\n",
      "The idea of this method is to maximize $E_q[ln(p(\\mathbf{f}))] + E_q[ln(p(\\mathbf{y}|\\mathbf{f}))] - E_q[ln(q(\\mathbf{f}|\\mathbf{y}))]$ with respect to $\\mu$ and $\\Sigma$, by reparametrizing $\\Sigma$ to reduce the dimension of variables to be optimized.\n",
      "However, such parametrization does not preserve convexity (concave in this setting) according to <a href=\"http://arxiv.org/abs/1306.1052\">the paper</a> of Khan et. al. in 2013. \n",
      "Except for this method, convexity holds in other variational methods implemented in Shogun. We demonstrate how to apply this variational method in Shogun to the USPS example."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "likelihood = LogitVGLikelihood()\n",
      "#using the default linesearch method\n",
      "learning_example2(KLCovarianceInferenceMethod, likelihood, x_train, x_test, y_train, y_test)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remark: This method may be the first variational method used in GPC models. However, in practice, it is the slowest method among all implemented variational methods in Shogun."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Mean-field Variational Method in GPC Models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This method is also called the factorial variational method in <a href=\"http://jmlr.org/papers/volume9/nickisch08a/nickisch08a.pdf\">the paper</a> of Nickisch and Rasmussen in 2008.\n",
      "The idea of mean-field variatonal method is to enforce artificial structure on the co-variance matrix, $\\Sigma$ of $q(\\mathbf{f})$.\n",
      "In mean-field variatonal method, $\\Sigma$ is a diagonal positive-definite matrix.\n",
      "We demonstrate how to apply this variational method in Shogun to the USPS example."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "likelihood = LogitVGLikelihood()\n",
      "#using the default linesearch method\n",
      "learning_example2(KLApproxDiagonalInferenceMethod, likelihood, x_train, x_test, y_train, y_test)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remark: This method could be as fast as the Laplace method in Shogun in practice but ignores all off-diagonal elements of the covariance matrix. However, in some case (ie, there are a lot of noise in data set), such structure regularization may bring some promising result compared to other variational methods, which can be observed in the previous session for the sonar data set."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Cholesky Variational Method in GPC Models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This method is proposed in <a href=\"http://jmlr.org/proceedings/papers/v15/challis11a/challis11a.pdf\">the paper</a> of Challis and Barber in 2011.\n",
      "The idea of this method is to maximize $E_q[\\ln(p(\\mathbf{f}))] + E_q[\\ln(p(\\mathbf{y}|\\mathbf{f}))] - E_q[\\ln(q(\\mathbf{f}|\\mathbf{y}))]$ in terms of the Cholesky representation, C, for the covariance matrix of $q(\\mathbf{f})$, where $\\Sigma=CC^T$ and C is a lower triangular matrix.\n",
      "We demonstrate how to apply this variational method in Shogun to the USPS data set as below."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "likelihood = LogitVGLikelihood()\n",
      "#using the default linesearch method\n",
      "learning_example2(KLCholeskyInferenceMethod, likelihood, x_train, x_test, y_train, y_test)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remark: This method may be faster to learn the complete structure of the covariance matrix from data than covariance variational method.\n",
      "The reason is solving linear system related to $\\Sigma$ is required at each optimization step.\n",
      "In Cholesky representation, the time complexity at each optimization step is reduced to O(n^2) from O(n^3) compared to covariance variational method.\n",
      "However, the overall time complexity is still O(n^3)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Dual Variational Method in GPC Models"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This method is proposed in <a href=\"http://arxiv.org/abs/1306.1052\">the paper</a> of Khan et. al. in 2013.\n",
      "The idea of this method is to optimize $E_q[ln(p(\\mathbf{f}))] + E_q[ln(p(\\mathbf{y}|\\mathbf{f}))] - E_q[ln(q(\\mathbf{f}|\\mathbf{y}))]$ in <a href=\"http://en.wikipedia.org/wiki/Duality_%28optimization%29\">Lagrange dual</a> form instead of the primal form used in covariance variational method via explicitly expressing constraint in the covariance matrix, $\\Sigma$, in terms of auxiliary variable.\n",
      "However, the time complexity of each optimization step is still O(n^3) and variational bound usually is required to approximate $E_{q_i}[ln(p(\\mathbf{y_i}|\\mathbf{f_i})]$ for GPC in this method.\n",
      "We demonstrate how to apply this variational method in Shogun to the USPS data set as below."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "likelihood = LogitDVGLikelihood()\n",
      "likelihood.set_strict_scale(0.1)\n",
      "#using the default linesearch method\n",
      "learning_example2(KLDualInferenceMethod, likelihood, x_train, x_test, y_train, y_test)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remark: This method requires O(N) variables to be optimized while the Cholesky method requires O(N^2) variables to be optimized. This fact is important for large-scale cases. A promising observation in practice is that this method is faster than covariance variational method and may not be slower than the Cholesky method. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Variational Bounds"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we discuss how to approximate $E_{q_i}[ln(p(\\mathbf{y_i}|\\mathbf{f_i})]$.\n",
      "Since this term is about one-dimensional integration, for GPC we can use <a href=\"http://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature\">Gauss\u2013Hermite quadrature</a> to appromxiate this term. In Shogun, the corresponding implementation of Bernoulli-logistic likelihood for variational inference is called <a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLogitVGLikelihood.html\">LogitVGLikelihood</a> . This likelihood can be used for all variational methods except the dual variational method.\n",
      "\n",
      "The piecewise variational bound proposed at <a href=\"http://www.icml-2011.org/papers/376_icmlpaper.pdf\">the paper</a> of Benjamin M. Marlin et. al. in 2011 is also implemented in Shogun, which is called <a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLogitVGPiecewiseBoundLikelihood.html\">LogitVGPiecewiseBoundLikelihood</a>. This likelihood can be used for all variational methods except the dual variational method.\n",
      "\n",
      "For dual variational method, we use the variational bound in <a href=\"http://www.cs.cmu.edu/~lafferty/pub/ctm.pdf\">the paper</a> of DM Blei et. al. in 2007. The corresponding implementation in Shogun is called <a href=\"http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CLogitVGPiecewiseBoundLikelihood.html\">LogitVGPiecewiseBoundLikelihood</a>. Note that when using LogitDVGLikelihood, the bound is only enabled for the dual variational method. When other variational methods use this class, it uses one-dimensional integration instead of the variational bound. For detailed discussion about variational bounds in GPC models, please refer to <a href=\"https://circle.ubc.ca/handle/2429/43640\">Khan's work</a>."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "likelihood = LogitVGPiecewiseBoundLikelihood()\n",
      "likelihood.set_default_variational_bound()\n",
      "likelihood.set_noise_factor(1e-15)\n",
      "inference_methods=[\n",
      "                  KLCholeskyInferenceMethod,\n",
      "                  KLApproxDiagonalInferenceMethod,\n",
      "                  KLCovarianceInferenceMethod,\n",
      "                  ]\n",
      "for inference in inference_methods:\n",
      "    #using the default linesearch method\n",
      "    learning_example2(inference, likelihood, x_train, x_test, y_train, y_test)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Optimizing Hyper-parameters"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We demonstrate how to optimize hypyer-parameters in Shogun for applying the mean-field variational method to the USPS data set."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def learning_example3(inference, likelihood, x_train, x_test, y_train, y_test):\n",
      "    \"\"\"\n",
      "    This function is used to train a GP classifer given an inference method and hyper-parameters are automatically optimized\n",
      "        \n",
      "    Parameters:\n",
      "        inference - an inference methods used to train the classifer\n",
      "        likelihood - likelihood function to model data\n",
      "        x_train, x_test - X for training and testing\n",
      "        y_train, y_test - Y for training and testing\n",
      "        \n",
      "    Returns:\n",
      "        Nothing\n",
      "    \"\"\"\n",
      "    error_eval = ErrorRateMeasure()\n",
      "    mean_func = ZeroMean()\n",
      "    kernel_log_sigma = 1.0\n",
      "    kernel_sigma = 2*exp(2*kernel_log_sigma);\n",
      "    kernel_func = GaussianKernel(10, kernel_sigma)\n",
      "\n",
      "    #sample by 1\n",
      "    labels_train = BinaryLabels(y_train)\n",
      "    labels_test = BinaryLabels(y_test)\n",
      "    #feature by sample\n",
      "    features_train=RealFeatures(x_train)\n",
      "    features_test=RealFeatures(x_test)\n",
      "\n",
      "    kernel_log_scale = 1.0\n",
      "\n",
      "    inf = inference(kernel_func, features_train, mean_func, labels_train, likelihood)\n",
      "    print \"\\nusing %s\"%inf.get_name()\n",
      "    \n",
      "    inf.set_scale(exp(kernel_log_scale))\n",
      "    try:\n",
      "            inf.set_lbfgs_parameters(100,80,3,80)\n",
      "    except:\n",
      "            pass\n",
      "\n",
      "    gp = GaussianProcessClassification(inf)\n",
      "\n",
      "    # evaluate our inference method for its derivatives\n",
      "    grad = GradientEvaluation(gp, features_train, labels_train, GradientCriterion(), False)\n",
      "    grad.set_function(inf)\n",
      "\n",
      "    # handles all of the above structures in memory\n",
      "    grad_search = GradientModelSelection(grad)\n",
      "\n",
      "    # search for best parameters and store them\n",
      "    best_combination = grad_search.select_model()\n",
      "\n",
      "    # apply best parameters to GP\n",
      "    best_combination.apply_to_machine(gp)\n",
      "\n",
      "    # we have to \"cast\" objects to the specific kernel interface we used (soon to be easier)\n",
      "    best_width=GaussianKernel.obtain_from_generic(inf.get_kernel()).get_width()\n",
      "    best_scale=inf.get_scale()\n",
      "    print \"Selected kernel bandwidth:\", best_width\n",
      "    print \"Selected kernel scale:\", best_scale\n",
      "\n",
      "    start = time.time()\n",
      "    gp.train()\n",
      "    end = time.time()\n",
      "    print \"cost %s seconds at training\"%(end-start)\n",
      "    nlz=inf.get_negative_log_marginal_likelihood()\n",
      "    print \"the negative_log_marginal_likelihood is %.4f\"%nlz\n",
      "    start = time.time()\n",
      "    #classification on train_data\n",
      "    pred_labels_train = gp.apply_binary(features_train)\n",
      "    #classification on test_data\n",
      "    pred_labels_test = gp.apply_binary(features_test)\n",
      "    end = time.time()    \n",
      "    print \"cost %s seconds at prediction\"%(end-start)\n",
      "    \n",
      "    error_train = error_eval.evaluate(pred_labels_train, labels_train)\n",
      "    error_test = error_eval.evaluate(pred_labels_test, labels_test)\n",
      "    \n",
      "    print \"Train error : %.2f Test error: %.2f\\n\" % (error_train, error_test);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "likelihood = LogitVGLikelihood()\n",
      "likelihood.set_noise_factor(1e-15)\n",
      "#using the default linesearch method\n",
      "learning_example3(KLApproxDiagonalInferenceMethod, likelihood, x_train, x_test, y_train, y_test)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Reference"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* Rasmussen, Carl Edward. \"Gaussian processes for machine learning.\" (2006).\n",
      "* Bishop, Christopher M. Pattern recognition and machine learning. Vol. 1. New York: springer, 2006.\n",
      "* Nickisch, Hannes, and Carl Edward Rasmussen. \"Approximations for binary Gaussian process classification.\" Journal of Machine Learning Research 9 (2008): 2035-2078.\n",
      "* Challis, Edward, and David Barber. \"Concave Gaussian variational approximations for inference in large-scale Bayesian linear models.\" International conference on Artificial Intelligence and Statistics. 2011.\n",
      "* Emtiyaz Khan, Mohammad, et al. \"Fast dual variational inference for non-conjugate latent gaussian models.\" Proceedings of The 30th International Conference on Machine Learning. 2013.\n",
      "* Marlin, Benjamin M., Mohammad Emtiyaz Khan, and Kevin P. Murphy. \"Piecewise Bounds for Estimating Bernoulli-Logistic Latent Gaussian Models.\" ICML. 2011.\n",
      "* Khan, Mohammad. \"Variational learning for latent Gaussian model of discrete data.\" (2012).\n",
      "* Wang, Chong, and David M. Blei. \"Variational inference in nonconjugate models.\" The Journal of Machine Learning Research 14.1 (2013): 1005-1031.\n",
      "* Paisley, John, David Blei, and Michael Jordan. \"Variational Bayesian inference with stochastic search.\" ICML. 2012."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Soon to Come"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* Large-scale sparse variational inference in GPC models\n",
      "* Stochastic variational inference in GPC models"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}